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Transitioning from MODIS to VIIRS: an analysis of inter-consistency 


of NDVI data sets for agricultural monitoring 


The Visible/Infrared Imager/Radiometer Suite (VIIRS) aboard the Suomi 
National Polar-orbiting Partnership (S-NPP) satellite was launched in 2011, in 
part to provide continuity with the Moderate Resolution Imaging 
Spectroradiometer (MODIS) instrument aboard National Aeronautics and Space 
Administration’s (NASA) Terra and Aqua remote sensing satellites. The VIIRS 
will eventually replace MODIS for both land science and applications and add to 
the coarse-resolution, long term data record. It is, therefore, important to provide 
the user community with an assessment of the consistency of equivalent products 
from the two sensors. For this study, we do this in the context of example 
agricultural monitoring applications. Surface reflectance that is routinely 
delivered within the M{O,Y}D09 and VNPO9 series of products provide critical 
input for generating downstream products. Given the range of applications 
utilizing the normalized difference vegetation index (NDVI) generated from 
M{O,Y}D09 and VNP09 products and the inherent differences between MODIS 
and VIIRS sensors in calibration, spatial sampling, and spectral bands, the main 
objective of this study is to quantify uncertainties related the transitioning from 
using MODIS to VIIRS-based NDVI’s. In particular, we compare NDVI’s 
derived from two sets of Level 3 MYD09 and VNP09 products with various 
spatial-temporal characteristics, namely 8-day composites at 500 m spatial 
resolution and daily Climate Modelling Grid (CMG) images at 0.05° spatial 
resolution. Spectral adjustment of VIIRS I1 (red) and I2 (near infra-red — NIR) 
bands to match MODIS/Aqua b1 (red) and b2 (NIR) bands is performed to 
remove a bias between MODIS and VIIRS-based red, NIR, and NDVI estimates. 
Overall, red reflectance, NIR reflectance, NDVI uncertainties were 0.014, 0.029 
and 0.056 respectively for the 500 m product and 0.013, 0.016 and 0.032 for the 
0.05° product. The study shows that MODIS and VIIRS NDVI data can be used 
interchangeably for applications with an uncertainty of less than 0.02 to 0.05, 
depending on the scale of spatial aggregation, which is typically the uncertainty 


of the individual dataset. 
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1. Introduction 


The Moderate Resolution Imaging Spectroradiometer (MODIS) aboard the Terra and 
Aqua remote sensing satellites has been successfully imaging the Earth’s surface since 
2000 and 2002, respectively. With spatial resolutions of 250 m, 500 m and 1 km and 36 
spectral bands, MODIS provides temporal composites and daily images with near real- 
time access to data (Davies et al. 2015) that are critical to many applications. The 
portfolio of the MODIS-based land products has been expanding and improving through 
5 Collections, to include surface reflectance (SR), vegetation indices (VIs), biophysical 
parameters (leaf area index and fraction of absorbed photosynthetically active 
radiation), net and gross primary productivity, bidirectional reflectance distribution 
function (BRDF), albedo, temperature and land cover. 

Surface reflectance that is delivered within the M{O,Y}D09 series of products 
(Vermote and Kotchenova 2008) provides a critical input for generating such 
downstream products and needs to be of the highest possible quality, so that minimal 
uncertainties propagate in the dependent/downstream products. An important 
downstream product is the normalized difference vegetation index (NDVI). NDVI has 
been one of the most important and widely applicable VI’s dating back to the Advanced 
Very High Resolution Radiometer (AVHRR) instruments aboard National Oceanic and 
Atmospheric Administration (NOAA) satellites (Justice et al. 1985), and is used in 
various agricultural applications including crop yield prediction (Becker-Reshef et al. 
2010a; Franch et al. 2015, 2017; Johnson 2016; Kogan et al. 2014; Meroni et al. 2016), 
crop mapping (Chang e¢ al. 2007; Pittman et al. 2010; Skakun et al. 2017; Xiao et al. 
2005), crop calendar and phenology analysis (Sakamoto et al. 2010; Whitcraft et al. 
2015), and drought monitoring and crop state assessment (AghaKouchak et al. 2015; 


Gu et al. 2007; Karl et al. 2012). More importantly, VIs derived from MODIS surface 


reflectance products have been integrated into operational agricultural monitoring 
systems at global, national and regional scales (Becker-Reshef et al. 2010b). With the 
MODIS Terra sensor already experiencing degradation (Wang et al. 2012), it is 
advisable to establish continuity observations for these applications. 

The Visible/Infrared Imager/Radiometer Suite (VIIRS) aboard the Suomi 
National Polar-orbiting Partnership (S-NPP) satellite was launched in 2011 and was 
planned to provide continuity with MODIS (Justice et al. 2013). VIIRS images the 
Earth’s surface in 22 spectral bands at 375 m (I bands) and 750 m (M bands) spatial 
resolution. A series of VIIRS-based surface reflectance products VNPO9, analogous to 
the M{O,Y}D09 suite, is being routinely generated (Vermote et a/. 2014) using the 
same approach for atmospheric correction as for MODIS (Vermote and Kotchenova 
2008; Vermote et al. 2002). VIIRS will eventually replace MODIS for both land science 
and applications and add to the coarse-resolution, long term data record. It is therefore, 
important to provide the user community with an assessment of the consistency of 
equivalent products from the two sensors. For this study, we do this in the context of 
example agricultural monitoring applications. Previous studies have provided some 
insight into continuity and inter-comparison issues between MODIS and VURS using 
simulated data (Fan and Liu 2016; Kim et al. 2010; Van Leeuwen et al. 2006; Miura et 
al. 2013), top-of-atmosphere (TOA) NDVI and top-of-canopy (TOC) enhanced 
vegetation index (EVI) (Fan and Liu 2017; Obata et al. 2016; Vargas et al. 2013) and 
AERONET-based validation (Shabanov et al. 2015). Given the range of applications 
utilizing NDVI generated from M{O,Y}D09 and VNP09 products and the inherent 
differences between MODIS and VIIRS instruments in terms of calibration, spatial 
sampling and spectral bands, the main objective of this study is to quantify uncertainties 


related to transitioning from MODIS to VURS based NDVI’s. In particular, we compare 


NDVI derived from two sets of Level 3 products (MYD09 and VNPO9) with different 
spatial-temporal characteristics. For this study, we selected: (i) 8-day composited 
products at 500 m spatial resolution, as composited data are commonly used in 
agricultural applications to minimize the impact of cloud cover and (ii) daily Climate 
Modelling Grid (CMG) images at 0.05° spatial resolution, as increasingly, daily data are 
being used to avoid losing high temporal frequency good observations eliminated by 
temporal compositing (Franch et al. 2017). The comparison is performed particularly 
for the MYDO09 products, with similar afternoon overpass times from Aqua and S-NPP 


satellites (i.e. 13:30 local time). 


2. Surface Reflectance Products M{O,Y}D09 and VNP09 


The M{O,Y}D09 (Vermote et al. 2015) and VNPO9 (Roger et al. 2016) products suites 
provide an estimate of the surface spectral reflectance for the corresponding MODIS 
and VIIRS spectral bands, as would have been measured at ground level if there were 
no atmospheric scattering or absorption. The same atmospheric correction algorithm 
which uses the Second Simulation of a Satellite Signal in the Solar Spectrum, Vector 
(6SV) radiative transfer code and internal algorithm for aerosol retrieval is applied to 
both MODIS and VIIRS (Vermote et al. 2014; Vermote and Kotchenova 2008). 
Corrections are made for the effects of molecular gases, including ozone and water 
vapour, and for the effects of atmospheric aerosols. 

M{O,Y}D09 is a seven-band product computed from the MODIS Level 1B 
bands | to 7. VNPO9 is a twelve-band product computed from the Land SIPS V1 Level 
1B bands I1-I3, M1—M5, M7, M8, M10, and M11. Both M{O,Y}D09 and VNP09 
include daily Level 2G (L2G) data, that have been mapped to the sinusoidal grid, and 


Level 3 (L3) data, that have been spatially and/or temporally aggregated. Table 1 and 2 


provide details on Level 2G and Level 3 products from the M{O,Y}D09 and VNP09 


series. 


[Table 1 near here] 


[Table 2 near here] 


For the temporal compositing process, each pixel containing the single best 
possible L2G observation during an 8-day period (hereafter referred as ‘best pixel’) is 
selected on the basis of high observation coverage, low sensor angle, the absence of 
clouds or cloud shadow, and aerosol loading. For spatial aggregation to CMG grid, an 
area weighted average of the best quality observations from the L2G product is used. 
The CMG product also provides the number of 250 m (for bands 1-2) and 500 m (for 


bands 3-7) best quality pixels which were used for averaging at 0.05° spatial resolution. 


3. Methodology 


NDVI products from MODIS and VIIRS at different spatial and temporal resolution 
were compared in this study. In particular, the comparison was performed on a per-pixel 
basis for MYD09 and VNPO9 products at 500 m and 0.05° (CMG) resolution, 
respectively. No aggregation (within a window) was performed for either resolution as 
the goal was to compare products at their ‘native’ resolutions. Bands b1 (red) and b2 
(near infra-red (NIR)) from MODIS and I1 (red) and I2 (NIR) from VIIRS were used to 
calculate the NDVI using a standard formula (Tucker 1979): (NDVI) = (pnir- 
Pred)/(PNIR+Pred), Where pnir and Pred are Surface reflectance values in NIR and red 
spectral bands, respectively. The VIIRS I1 and I2 bands were used instead of M5 and 
M7 bands for the CMG product because the spectral response functions from the I 


bands are more similar to those from MODIS, especially in the red (Figure 1). 


[Figure 1 near here] 


The 500 m 8-day composite product, comparison was performed for four tiles of 
the MODIS sinusoidal grid (h10v04, h10v05, h11v04 and h11v05), covering the 
Midwest US (Corn Belt, Figure 2), which is a major agricultural production region in 


the U.S. 


[Figure 2 near here] 


Though MODIS/Aqua and VIIRS sensors image the Earth’s surface at 
approximately the same time of day, the day of the year (DOY) within an 8-day period, 
for which the ‘best pixel’ value is selected for MYDO9A land VNPO9A1 products, 
might be different. Therefore, only the same DOY observations were used for 
comparison. In addition, only close to nadir observations from both sensors, 1.e. with 
view zenith angle (VZA) less than 7.5°, were considered to reduce the effects of spatial 
resolution and BRDF. 

A comparison of CMG products, namely MYDO9CMG and VNPO9CMG, was 
performed globally for land pixels. Daily MODIS and VIIRS CMG products exhibit 
different viewing geometries, and therefore BRDF correction is necessary to normalize 
the surface reflectance values. For this, we applied the VJB algorithm (Vermote et al. 
2009) for both MODIS/Aqua and VURS. The surface reflectance values in red and NIR 
bands from MODIS and VIIRS were normalized to the (45°, 0°) solar and viewing 


angles: 


1+VF, (45,0,0)+RF> (45,0,0) 
14+VF,(65,0y,9)+RF2(05,0y,0)’ 


p\(45,0,0) = p(Os, y, P) (1) 


where 9, is the solar zenith angle; 0, is the sensor VZA; ¢ is the relative azimuth angle; 
F, is the volume scattering kernel, based on the Rossthick function, but corrected for the 
Hot-Spot process; Fz is the geometric kernel, based on the Li-sparse reciprocal function; 
V and R are free parameters that are estimated for each pixel at CMG resolution using 
the BRDF inversion technique. We refer the reader to (Vermote et al. 2009) for the 
details of the VJB algorithm implementation. 

Since MODIS and VIIRS spectral response functions in red and NIR bands 
exhibit differences (Figure 1), corresponding spectral adjustments should be performed 
to reduce differences in red, NIR and NDVI estimates derived from MODIS and VIIRS 
sensors. In this study, surface reflectance values from VIIRS were adjusted to those of 
MODIS/Aqua. For this, corresponding relationships between red and NIR bands from 


MODIS and VURS were developed using the following equations: 


Pred = Gap ae bredPNIR> (2) 


PNIR = QninPrea + DyirPNiR: (3) 


where pia, PNIR> Prea> Pir ate surface reflectance values in red and NIR for MODIS 
(superscript M) and VIIRS (superscript V), and Qyeq, Drea, Qnir> DniR are conversion 
coefficients estimated from data using the ordinary least squares (OLS) regression. Note 
that these relationships are without the constant term to ensure that ‘black’ surfaces 
have the same reflectance values for both sensors. It is also expected that both the sums 
of coefficients, namely Qpeg + Dreqg and AyiR + Dyjip, will be close to | to ensure 
continuity of reflectance values for ‘bright’ surfaces such as clouds. Spectral adjustment 
for red and NIR bands for MYDO9A1/VNPO9A1 (500 m) and 
MYDO9CMG/VNPO9CMG (0.05°) products using Equations (2)-(3) was performed on 


a yearly basis and for the entire period of 2012-2016, to analyse the ‘temporal stability’ 


of the conversion coefficients. 

The NDVI anomaly is another indicator often used to analyse how the current 
vegetation condition relates to that of the previous years (AghaKouchak et al. 2015; 
Becker-Reshef et al. 2010b; Gu et al. 2007; Karl et al. 2012; Meroni et al. 2016). Here, 
we calculate a multi-year median of NDVI for each pixel from MODIS/Aqua, and 
combined MODIS/Aqua and adjusted VURS at CMG resolution. More specifically, for 
MODIS/Aqua data sets only, a median NDVI is calculated for 2002-2012. For a 
combination of MODIS and VIIRS NDVI data, a median NDVI is calculated for a set 
of NDVI values concatenated from MODIS/Aqua (2002-2011) and adjusted VIIRS 
(2012-2016). Therefore, we compare two cases: when VIIRS adjusted data, starting 
from 2012, are used to update the median NDVI values from MODIS/Aqua, and when 


only MODIS/Aqua is used to calculate the median NDVI (2002-2016): 
yo dian = Median({yM|t = 2002..2016}), (4) 
ymedian = Median({y|t = 2002. .2011}, {yf |t = 2012..2016}), (5) 


where ¢ is the time, expressed here in years, y{’ and y/ are MYDO9CMG and 
VNPO9CMG derived NDVI values for MODIS/Aqua and VIIRS (adjusted with 
Equations (2)-(3)), respectively. 

To quantify the differences and uncertainties between MYDO9 and VNP09 
products and NDVI-derived estimates, a standard APU analysis is performed (Vermote 


and Kotchenova 2008) with the following set of metrics: 


* accuracy (A) that shows the average bias between estimates 


A= “YN (7 - yl") (6) 


* precision (P) that shows repeatability of the estimates 


1 2 
P= ar = 9 =A) (7) 
* uncertainty (U) that is the root mean square error 
ee V _.M)2 8 
aeey ire Die1(y7 — yi") (8) 


where yj’ and yM are VIIRS and MODIS derived values (surface reflectance or NDVI), 


respectively, and N is the number of values to be compared. 


4. Results 


4.1. Comparison of MYDO9AI and VNPOQAI derived NDVI at 500 m resolution 


Figure 3 shows the difference between the DOY selected within the 8-day period for 
MYDO9A1 and VNPO9A1 products for all land pixels (excludes water pixels) in tiles 
h10v04, h10v05, h11v04 and h11v05 for 2012-2016. The same DOY is selected in 
35.4% cases; in 55.6% cases, the DOY difference is | day or less; and in 31% cases, the 
DOY difference is 3 to 7 days. Therefore, when utilizing these products jointly, one has 
to take into account differences in surface reflectance values or derived VI’s caused by 


the DOY selected and possible surface changes within the compositing period. 


[Figure 3 near here] 


The distribution of MODIS/Aqua VZAs in MYDO9A1 and difference between 


VZAs for MYDO9A1 and VNPO9A1 products for the four tiles in US and the 2012- 


2016 period, are shown in Figure 4 and Figure 5, respectively. For MODIS/Aqua VZAs 


distribution, in 27.3% cases the VZA is 0° to 10°; in 35% of cases, VZA is more than 
30° which translates to the effective along-scan spatial resolution of more than 850 m 
(Figure 4). In 38% and 82% of cases, the difference between MODIS/Aqua and VIIRS 
VZAs within the 8-day composites is within -10° to 10° and -30° to 30°, respectively; in 


18%, the absolute difference is more than 30°. 


[Figure 4 near here] 


[Figure 5 near here] 


The results of ‘temporal’ stability of conversion coefficients are presented in 


Table 3. 


[Table 3 near here] 


In general, there is a temporal ‘stability’ for coefficients a,¢q and by;pr as the 
coefficient of variation (CV) is around 1% which is consistent with the calibration 
performance (Vermote et al. 2014). Coefficients b-eg and Ayjp provide a maximum 2% 
and 6% contribution to the surface reflectance values for red and NIR bands. Less 
stability is observed for these coefficients, b-eq and Ayr, with CV of 6.8% and 44.3%, 
respectively. It should be noted that larger deviations are observed for the initial years 
of VIIRS operation, namely 2012 and 2013, while relatively better stability is observed 
for 2014-2016. 

The derived coefficients for 2012—2016 (Table 3) were used to adjust red and 


NIR reflectance values from VHURS to match the MODIS ones, and compute the NDVI. 


Comparisons of the VIIRS and MODIS/Aqua derived NDVI values with and without 


spectral adjustment are shown in Figure 6 and Figure 7, respectively. 


[Figure 6 near here] 


[Figure 7 near here] 


Spectral adjustment removed the bias between MODIS/Aqua and VIIRS-derived 
red, NIR and NDVI values (Figure 6 and Figure 7). Overall red reflectance, NIR 
reflectance and NDVI uncertainties were 0.014, 0.029 and 0.056, respectively, when 
considering the same day observation pixels (Figure 6). These uncertainties increased to 
0.018, 0.034 and 0.064 (a 14% increase), respectively, when the absolute difference 
between DOY for ‘best pixels’ from MODO9A1 and VNPO9A1 products was 3 to 7 


days while VZAs for both sensors were less than 7.5° (Figure 8). 


[Figure 8 near here] 


Uncertainties can be further reduced when NDVI values at 500 m resolution are 
spatially aggregated. For example, within the agriculture application domain, NDVI is 
usually averaged over administrative regions to correlate with crop yield values 
(Becker-Reshef et al. 2010a; Franch et al. 2015, 2017; Johnson et al. 2016; Kogan et al. 
2014). Figure 9 shows an example of such aggregation: NDVI values derived from 
MODO9A1 and VNPO9A1 products were averaged for Harper County in Kansas (US) 
for 500 m pixels with a winter wheat proportion larger than 50%. Winter wheat 
proportions were derived from USDA CDL maps for 2012—2016 (Johnson and Mueller, 


2010). The spatial aggregation decreased uncertainties to 0.021 (2.67 times), compared 


to a per-pixel (at 500 m) derived uncertainty of 0.056. Spectral adjustment reduced the 


bias (accuracy) from 0.017 to -0.003. 


[Figure 9 near here] 


4.2. Comparison of MYDO9CMG and VNP09CMG derived NDVI at 0.05° 


resolution 


Table 4 shows the derived coefficients from Equations (2)—(3) of the regressions to 
adjust VIIRS red (11) and NIR (12) surface reflectance values to MODIS using yearly 
data and the whole 2012—2016 period. Compared to the 500 m products (MYDO9A1 


and VNP09A1), better temporal ‘stability’ is observed at CMG resolution. 


[Table 4 near here] 


The conversion coefficients from Table 4 derived for 2012—2016 were used to 
adjust VIIRS I1 (red) and I2 (NIR) bands and to compute NDVI. Figure 10 shows 
comparison of daily NDVI values at 0.05° resolution for all land pixels for 2012—2016 


(almost 2 x 10° CMG pixels). 


[Figure 10 near here] 


The spectral adjustment removed the bias between MODIS/Aqua and VIURS 
derived NDVI, and the resulting uncertainty for NDVI was 0.032, 1.75 times smaller 
than for the 500 m products. Uncertainties for red and NIR spectral bands were 0.013 


and 0.016, respectively. CMG based data were spatially aggregated over Harper County 


to provide a daily NDVI time series for winter wheat (Figure 11). We selected the top 
5% purest winter wheat pixels at CMG resolution to calculate NDVI, the same way it 
was done for the generalized empirical winter wheat yield forecasting model (Becker- 


Reshef et al. 2010a). 


[Figure 11 near here] 


With aggregation at Harper County scale, both 500 m and CMG products 
yielded similar uncertainty of 0.022 when comparing MODIS/Aqua and VIIRS derived 
NDVI values (Figure 9 and Figure 11); however, the CMG derived NDVI time series is 


much denser thanks to daily observations. 


4.3. Comparison of MYDO9CMG and VNP09CMG derived NDVI anomalies at 


0.05° resolution 


Figure 12 shows a comparison of long-term median NDVI values calculated for 
MODIS/Aqua only (Equation (4)) and a combination of MODIS/Aqua and adjusted 
VIIRS (Equation (5)) for 2002—2016. Thanks to VIIRS adjustment, the bias is close to 
zero (-0.003), and the corresponding NDVI uncertainty is 0.03. Comparison of NDVI 
anomalies derived from BRDF corrected MOD09CMG and VNPO9CMG products is 


shown in Figure 13. The resulting uncertainty was found to be 0.033 at global scale. 


[Figure 12 near here] 


[Figure 13 near here] 


An example of NDVI values and medians for Iowa (US) is shown in Figure 14, 


and geographical distribution of NDVI anomalies from MODIS/Aqua and VIIRS for the 


same region is shown in Figure 15. Figure 15 shows good spatial consistency between 
and similar spatial patterns for NDVI anomalies computed from MODIS/Aqua and 


VIIRS sensors. 


[Figure 14 near here] 


[Figure 15 near here] 


5. Discussion 


NDVI is a widely used remote sensing based product which is used in several 
agricultural monitoring applications. Having high-quality long term NDVI data records 
is extremely important for studying spatiotemporal changes in Earth’s surface 
dynamics. This requires integration of data records from multiple sensors, including 
MODIS and VIIRS. VIIRS provides continuity to MODIS, and therefore it is important 
to enable a proper transition between products from these sensors, so the VIIRS based 
products can be ingested into existing MODIS based applications, and corresponding 
uncertainties are quantified and known. This study focused on comparing NDVI derived 
from the MYD09A1/VNPO9A1 and MYDO9CMG/VNPO9CMG surface reflectance 
products that provide a trade-off in terms of spatial (500 m versus 0.05°) and temporal 
resolution (8-day versus daily). In particular, MYDO9A1 and VNPO9A1 provide data at 
500 m resolution at the expense of temporal resolution (8-day). Although, through the 
compositing process, only high-quality pixels are selected, there are several differences 
between MODIS/Aqua and VIIRS-based products that influence the inter-consistency 
of the data sets. First and foremost, differences in spectral response functions in red and 
NIR bands of MODIS and VURS sensors (Figure 1) introduce a bias in surface 


reflectance values and NDVI estimations that can be removed through spectral 


adjustments. We found also, that only in 35.4% cases the DOY of the ‘best pixel’ within 
the 8-day period is the same for MYDO9A1 and VNPO9A1 products. Uncertainties of 
MYDO9A1 and VNPO9A1 derived NDVI values can be increased more than 14% when 
the difference between DOY increases to 7 days. Differences were also observed in 
VZAs for ‘best pixels’ in these products. The off-nadir VZAs values introduce two 
major issues: a reduced effective spatial resolution of MODIS with the increase of VZA 
(Figure 4) (thanks to the aggregation, this is not the case for VIRS (Campagnolo et al. 
2016; Pahlevan et al. 2017)), and BRDF effects for both MODIS and VIIRS. Therefore, 
users are encouraged to take these into consideration when developing applications at 
‘native’ 500 m resolution. Overall, the uncertainties between MYDO9A1 and VNPO9A1 
derived red reflectance, NIR reflectance and NDVI estimates at 500 m resolution for the 
same day and close to nadir (VZA<7.5°) observations for US for 2012—2016 were 
found to be 0.014, 0.029 and 0.056, respectively, with VIIRS to MODIS/Aqua spectral 
adjustment. 

A better consistency between MODIS/Aqua and VIIRS derived NDVI’s was 
observed at CMG scale at 0.05° resolution. Comparison of more than 2 x 10° global 
CMG pixels for 2012-2016 yielded red reflectance, NIR reflectance and NDVI 
uncertainties of 0.013, 0.016 and 0.032 respectively, after BRDF correction of 
MYDO9CMG and VNPO9CMG surface reflectance values with the VJB approach 
(Vermote et al. 2009), and spectral adjustment of VIIRS to MODIS/Aqua. 
Corresponding conversion coefficients for adjusting BRDF corrected VIIRS I1 (red) 
and I2 (NIR) bands to MODIS/Aqua b1 (red) and b2 (NIR) bands were calculated and 
showed good temporal stability within the 2012—2016 period. While the CMG based 
products provide a lower spatial resolution as compared to the 500 m products, they 


provide daily data that might be critical to applications, and, as it is shown in this study, 


better consistency between MODIS/Aqua and VIIRS with lower uncertainties. 
Uncertainties can be further reduced to 0.022 when NDVI values extracted from 500 m 
or CMG resolution products are spatially aggregated for administrative regions, as the 
derived NDVI can, for example, be correlated with crop yields (Becker-Reshef et al. 
2010a; Franch et al. 2015, 2017; Johnson et a/. 2016; Kogan et al. 2014). 

These results have certain implications when ingesting VIIRS data into existing 
MODIS-based models for agricultural monitoring, e.g. crop state assessment or crop 
yield modelling and forecasting. VIIRS data should be spectrally adjusted to match 
MODIS data, so no bias will propagate into the final estimates. Consider a model with 
crop yield linearly depending on the MODIS-based NDVI: y = a*(NDVD). Directly 
applying the VIIRS-based NDVI’s without spectral adjustment will result in higher crop 
yield estimates, since VIIRS-based NDVI is higher than the MODIS/Aqua-based NDVI 
(Figure 9 and Figure 11). For example, a slope between winter wheat yield (t ha”') and 
MODIS-based NDVI for Harper County in Kansas (US) was found to be 5.34 (Becker- 
Reshef et al. 2010a), while the VIIRS-based NDVI’s are on average 0.018 higher than 
the MODIS-based NDVI’s (Figure 11). Therefore, in such a case, the VIIRS-based 
winter wheat yield estimates will be on average 0.1 t ha’! higher than those from 
MODIS without spectral adjustment. 

Even with the bias removed, differences still exist between MODIS and VIIRS 
derived NDVI’s, and quantifying inter-consistency between the sensors can be helpful 
in providing the final error of estimates for NDVI-based agricultural products. Consider 
again the example of Kansas (US) where the winter wheat yield model is estimated to 
have an RMSE error of 0.18 t ha’! or 7% (Becker-Reshef et al. 2010a). When applying 
the VIIRS-based NDVI’s to the model, due to the MODIS-VIIRS NDVI uncertainty of 


0.022 (Figure 11), the error of winter wheat yield estimates will be 


0.182 + (5.34x0.022)2 = 0.22 tha’ or 8.5%. Therefore, inconsistencies between 
VIIRS and MODIS based NDVI’s will lead to the increase of resulting crop yield 
uncertainties. 

In terms of NDVI anomalies, median values are calculated from a sufficiently 
long data record to identify ‘normal’ vegetation conditions, but not so long that the land 
use or cropping system being observed have changed significantly. Figure 16 shows the 
timeline of the three remote sensing satellites imaging the Earth’s surface with an 
afternoon overpass, that will be used to form the land long term record: MODIS/Aqua, 
S-NPP and Joint Polar Satellite System (JPSS). It is expected that MODIS/Aqua will 
continue its nominal operations until 2022 (personal communication, Robert Wolfe, 
NASA Goddard Space Flight Center, June 2017) and JPSS-1 is planned to be launched 
at the end of 2017. At the time of writing, MODIS/Aqua has a 15 year data record 
which is used to calculate the median NDVI value. At the MODIS/Aqua end of life 
(2022), the data record would be 20 years and the VIIRS/S-NPP record would be 10 
years. Inter-use of data products from these sensors is therefore likely to continue to be 
desirable; however, if the NDVI data records are combined, one should do so with an 


awareness of NDVI anomaly inter-consistency uncertainties of 0.033. 


[Figure 16 near here] 


6. Conclusion 


The main focus of this study was to quantify uncertainties between MODIS/Aqua and 
VIIRS based NDVI calculated from the suite of MYD09 and VNPO9 products that 


provide an estimate of surface reflectance, which provides the basis for the NDVI. 


Because of differences in spectral response functions in red and NIR bands of MODIS 
(b1, b2) and VURS (11, 12), there is a bias when comparing MODIS and VIURS 
estimates that is removed through a spectral adjustment procedure. Corresponding 
coefficients were calculated for MYD09 and VNPO9 products at 500 m (for US) and 
CMG (globally) spatial resolution using observations from 2012 to 2016 that can be 
further used by the user community in their research activities. At 500 m spatial 
resolution and 8-day temporal resolution, uncertainty between NDVI derived from 
MODIS/Aqua and VIIRS was 0.056 for the same day and close to nadir observations 
(VZA<7.5°). For daily BRDF corrected NDVI and NDVI anomalies values at 0.05° 
resolution, uncertainty was 0.032 and 0.033, respectively. Uncertainty between 
MODIS/Aqua VIIRS derived NDVI can be further reduced to 0.022 when aggregating 
NDVI values over administrative regions. The derived NDVI uncertainties for different 
MODIS/Aqua and VIIRS products can be used by user community to quantify 
uncertainties for high level products. With the launch of JPSS-1 VIIRS later this year, 
there will be a need for additional product inter-comparisons similar to this study, in the 


context of data inter-use and land long-term data records. 
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Table 1. The M{O,Y}D09 Collection 6 Product Suite (Vermote et al. 2015). 


Product name Product description 

(Terra / Aqua) 

MOD09GQ / MYD09GQ Surface Reflectance Daily L2G Global 250 m 
(bands 1, 2) 

MODO09GA / MYDO9GA Surface Reflectance Daily L2G Global 500 m and 1 km 
(bands 1-7) 

MOD09Q1 /MYDO9QI1 Surface Reflectance 8-Day L3 Global 250 m 
(bands 1, 2) 

MODO9A1 / MYDO9A1 Surface Reflectance 8-Day L3 Global 500 m 
(bands 1-7) 

MODO09CMG / MYDO9CMG | Surface Reflectance Daily L3 Global 0.05° CMG 
(bands 1-7) 


Table 2. The VNP09 Collection 1 Product Suite (Roger et a/. 2016). 


Product name Product description 


VNPO9GHKI Surface Reflectance Daily L2G Global 500 m 


(bands I1—I3) 


VNPO9IGIKI Surface Reflectance Daily L2G Global | km 


(bands M1—M5, M7, M8, M10, M11) 


VNPO9IGA Surface Reflectance Daily L2G Global 500 m and 1 km 


(bands I1-I3 (500m), M1—MS5, M7, M8, M10, M11 (1 km)) 


VNPO9A1 Surface Reflectance 8-Day L3 Global 500 m 


(bands I1—I3) 


VNPO9H1 Surface Reflectance 8-Day L3 Global 1 km 


(bands M1—M5, M7, M8, M10, M11) 


VNPO9CMG 


Surface Reflectance Daily L3 Global 0.05° CMG 


(bands I1—-I3, M1—M5, M7, M8, M10, M11) 


Table 3. Estimated conversion coefficients for spectral adjustment of red and NIR 
spectral bands from VNPO9A1 to MYDO9A1. Only pixels with the same DOY and 
close to nadir observations (VZA<7.5°) for MYDO9A1 and VNPO9A1 were considered. 


Period Qred Drea QNIR byir 
2012 0.9788 0.0174 0.0834 0.9394 
2013 0.9704 0.0185 0.0778 0.9417 
2014 0.9628 0.0204 0.0357 0.9622 
2015 0.9691 0.0196 0.0378 0.9622 
2016 0.9562 0.0176 0.0369 0.9533 
Mean + standard | 0.9674+0.0085 | 0.0187+0.0013 | 0.0543+0.0241 | 0.9518+0.0109 
deviation (0.9%) (6.8%) (44.3%) (1.1%) 
(coefficient of 

variation, %) 

2012-2016 0.9687 0.0184 0.0544 0.9518 


Table 4. Estimated coefficients for spectral adjustment of red and NIR spectral bands 


from MYDO9CMG and VNPO9CMG. 


Period Ared Dred QANnIR byir 
2012 0.9805 0.0187 0.0011 0.9720 
2013 0.9812 0.0190 0.0008 0.9733 
2014 0.9799 0.0189 0.0011 0.9724 
2015 0.9809 0.0184 0.0010 0.9730 
2016 0.9842 0.0143 0.0010 0.9730 
Mean + standard 0.9813+0.0017 | 0.0178 0.0010+0.0001 | 0.9727+0.0005 
deviation (0.2%) +0.0020 (10.8%) (0.1%) 
(coefficient of (11.2%) 

variation, %) 

2012-2016 0.9814 0.0178 0.0020 0.9717 
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Figure 1. Relative spectral response functions for MODIS/Aqua and VIIRS sensors in 
the red and NIR spectral domain. The functions for MODIS and VIIRS were derived 
from https://mcst.gsfc.nasa.gov/calibration/parameters and 


https://ncc.nesdis.noaa.gov/VIIRS/VURSSpectralResponseFunctions.php, respectively. 
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Figure 2. Illustration of four MODIS tiles over the US (in MODIS sinusoidal 
projection) used for comparison of 500 m 8-day composite products from MODIS 
(MYDO9A1) and VURS (VNP09A1) sensors. Shown also is a distribution of croplands 
derived from the USDA’s Cropland Data Layer (CDL) for 2016 (Johnson and Mueller 


2010). 
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Figure 3. Distribution of DOY difference for MYDO9A1 and VNPO9A1 8-day 
composite products at 500 m spatial resolution. All land pixels from MODIS tiles 
hl0v04, hl0v05, hl1v04 and hl1v05 for 2012-2016 were used to build the chart. 
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Figure 4. Distribution of VZA values for ‘best pixels’ selected within the 8-day period 
for the MYDO9A1 products. All land pixels from MODIS tiles h10v04, h10v05, h11v04 
and h11v05 for 2012—2016 were used to build the chart. Also shown are the effective 
MODIS along-scan and along-track pixel sizes for the nominal 500 m resolution 


depending on VZAs (Campagnolo and Montano 2014). 
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Figure 5. Distribution of the difference between VZA values for MYDO9A1 and 
VNPO9A1 8-day composite products at 500 m spatial resolution. All land pixels from 
MODIS tiles h10v04, hl0v05, hl1v04 and hl1v05 for 2012—2016 were used to build 
the chart. 
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Figure 6. A scatterplot of red, NIR and NDVI values derived from VNPO9A1 (after 
spectral adjustment) and MYDO9A1 at 500 m resolution (a), (c) (e). Corresponding 
APU analysis (5), (d) (f). Land pixels from MODIS tiles h10v04, h10v05, h11v04 and 
h11v05 for 2012—2016 and having the same DOY for MODIS and VIIRS and close to 
nadir observations (VZAs<7.5°) were considered. The light blues bars on (5), (d) (f) 
show the number of points used in each bin of surface reflectance or NDVI values from 
MODIS/Aqua (used as a reference). The APU values (Equations (6)-(8)) were 


computed for points in each bin and being shown in red (accuracy), green (precision) 


and blue (uncertainty). The pink represents the specified uncertainty based on 
theoretical error budget of the Collection 5 MODIS (Vermote and Kotchenova 2008): 


¥2(0.005 + 0.05p) for spectral bands and ¥2(0.02 + 0.02V/) for NDVI. The V2 term 


is used since we are focusing on inter-consistency of datasets and not validation. 
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Figure 7. The same as Figure 6 but without spectral adjustment. 
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Figure 8. The same as Figure 6 but the absolute difference between DOY for 
MYDO9A1 and VNPO9A1 is 3 to 7 days. 
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Figure 9. A time series of NDVI values derived from MYDO9A1 and VNPO9A1 8-day 


products at 500 m resolution for Harper County, one of the largest wheat producing 


counties in Kansas. Surface reflectance values in red and NIR bands, that were used to 


compute NDVI from VIIRS, were spectrally adjusted to match the MODIS/Aqua ones 


(using Equations (2)—(3) and derived coefficients from Table 3). Shown also are final 


winter wheat yields derived from USDA National Agricultural Statistics Service 


(NASS) statistics (a). The difference between aggregated NDVI values from 
MYDO9A1 and VNPO9A1 with and without spectral adjustment of the VIIRS bands 


(5). 
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Fig. 10. A scatterplot of red, NIR and NDVI values derived from VNPO9CMG (after 
spectral adjustment) and MYDO9CMG at 0.05° resolution (a), (c) (e). Corresponding 
APU analysis (5), (d) (f). The light blues bars on (6), (d), (f) show the number of points 
used in each bin of surface reflectance or NDVI values from MODIS (used as a 
reference). The APU values (Equations (6)-(8)) were computed for points in each bin 
and being shown in red (accuracy), green (precision) and blue (uncertainty). The pink 
represents the specified uncertainty based on theoretical error budget of the Collection 5 


MODIS NDVI (Vermote and Kotchenova 2008). 


0.9 + 
| 2.71 t/ha 2.07 t/ha 0.85 t/ha 1.86 t/ha 3.12 t/ha 
0.8 T x * 
| . t a 0 
07 3-4 * A bz 
| ‘ ° 7 . e e = 
S | + ‘ t . 1 & 
® 06 a iar ee eas, ane . i. op aaa ‘ i ie 7 gt 8 ges “* 
|p f ae ! r. | 
0.5 {e" “= etait | : s 1-}-- -- - Or} ft --- en a senate a= 
| . sal ‘ 8 if Hy! i " G . "Lee 7 . @ : ° 
| . s . . . - e . ‘ 
0.4 + ‘ S '] a P-|-~ 5 
| ; Aa f ti y j tee “f teen 
. . ' 
| 5 . 
03 + :ipt o.-- 24 * MYDO9CMG 
« VNPOSCMG 
0.2 + + } 
2012 2013 2014 Vear 2015 2016 (a) 
0.15 : 
” | 
3 | 
Z 0.10 +--------------------}--------- 22-22-2222 nn nnn nnn nnn nnn fen nnn nnn nnn nen epee enn n eee een e eee 
a} . 
5 * : 
: on ee oT ae | a 7 en aaa (ace oF ae A id eet a, ay s5--- 8-2 se 
c . . ot .* s we : ne we tee “5 “ SL J wt? eens 
6 * > ? *** * see een a) . x ag" 
“a A ee ee i re et A) oe [RP ee aS Bee | ge + +e 
« eztxe 2 8: 4 ie | * ee anr® cer. TAGS tcp! Beast amed + fa bod 
= 0.00 lew 5D ae ttt a eae Aa? HX ee pe hs ak 4 2 Pee a 
5 ao ot] Ps °. " + rn . oe ’ ohtag te ¢ ; 
g . e t ¢ of 
B -0.05 + --0---5--------- 2-2 penn nna nee nnn nn nn nn] on nn nn nn nnn nnn nen fen enn nnn nnn e eee ig Saeak ab aiateiaieiaieiaiaiaieiatee 
2 ot 
$ . 
© . ; 
0.10 | --4|----------- -- - : 
£ + With spectral adjustment (A=-0.006, P=0.021, U=0.022) 
a « Without spectral adjustment (A=0.018, P=0.021, U=0.028) 
-0.15 + A J 
2012 2013 2014 Year 2015 2016 (b) 


Figure 11. A time series of NDVI values derived from BRDF corrected MYDO9CMG 


and VNPO9CMG daily products at 0.05° resolution for Harper County. 


Surface 


reflectance values in red (11) and NIR (12) bands, that were used to compute NDVI 


from VURS, were spectrally adjusted to match the MODIS/Aqua ones. 


Shown also are 


final winter wheat yields derived from USDA NASS statistics (a). A difference between 
aggregated NDVI values from MYDO9CMG and VNPO9CMG with and without 


spectral adjustment for VURS bands (db). 
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Figure 12. Comparison of median NDVI values at CMG resolution derived from 
MODIS/Aqua (Equation (4)) and a combination of MODIS/Aqua and adjusted VIIRS 
(Equation (5)) for 2002-2016 at a daily timestamp. A combined (MODIS/Aqua and 
VIIRS) median NDVI versus MODIS/Aqua derived median NDVI is shown in (a); 


APU analysis (5). 
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Figure 13. Comparison of NDVI anomalies for 2016 at CMG resolution at daily 


timestamp. Anomalies were derived by subtracting daily NDVI values for 2016 from 
median values calculated for 2002—2015. Combined (MODIS/Aqua and VIIRS) NDVI 


anomalies versus MODIS/Aqua derived NDVI anomalies is shown in (a); APU analysis 


(b). 
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Figure 14. Corn growth dynamics derived from MODIS/Aqua and VURS in 2012 in 
Iowa (US) compared to the median NDVI values for 2002—2016 derived from 
MODIS/Aqua. Due to a drought, corn growth started to decrease significantly from 
June which resulted in a 25% yield reduction according to USDA NASS. 


93°W 92°~w S1°W 


NDVI Anomaly 
DE No vata 
| ery 
HB 040-03 
HB 030-02 
BR 0210-01 
DY 0.1 t0-0.025 
0.025 to 0.025 
PY) 0.025 100.1 
GB 0.1 002 
HM o2003 
HH oswo4 
HB osm06 


Figure 15. NDVI anomalies at 0.05° spatial resolution for the state of Iowa (US) derived 
from MODIS/Aqua (a), and adjusted VIIRS (b) data on 21 August 2012. Anomalies 
were computed by subtracting NDVI values from the median NDVI values for 2002— 
2016 derived from MODIS/Aqua. 
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Figure 16. Timeline of the three afternoon remote sensing satellites: Aqua, S-NPP and 


JPSS-1. 


